cap program drop frolichlogitdgp
program define frolichlogitdgp
syntax, design(integer)  n(integer) var(real) curves(numlist) [approx(string)]
qui{
	drop _all
	set obs `n'
	gen X=rnormal()
	summ X
	local Xmin=-5
	local Xmax=5
	gen x =(X-(`Xmin'))/((`Xmax')-(`Xmin')) //ranges from 0 to 1
	replace x=2*(x-0.5)                     //ranges from -1 to 1
	replace x=_pi*x                         //ranges from -pi to pi

	forvalues i=1(1)7 {
	  gen S`i'=sin(`i'*x)
	  gen C`i'=cos(`i'*x)
	}
	fouriercoeffs
	gen C0 = 1
	putmata X = (S1 S2 S3 S4 S5 C1 C2 C3 C4 C5 C0), replace
	mata: K = st_matrix("K")
    mata: Xk=X*K[.,`design']
	getmata Xk, replace
    gen u = runiform()
    replace u = log(u/(1-u))
    gen pX= exp(Xk)/(1+exp(Xk))
	gen Tstar= Xk + u
	gen T= Tstar>=0

	** Frolich Outcome Equation
    gen var = `var'
	gen e =(runiform()-.5)*sqrt(var)*sqrt(12)
	local X = "pX"
	foreach k of local curves{
		if `k'==1 gen m1 = 0.15 + 0.7* `X'                                     
		if `k'==2 gen m2 = 0.1 + `X'/2 + (1/2)*exp(-200*(`X'-0.7)^2)           
		if `k'==3 gen m3 = 0.8 - 2*(0.9-`X')^2 - 5*(`X'-0.7)^3 - 10*(`X'-0.6)^10 
		if `k'==4 gen m4 = 0.2 + sqrt(1-`X') - 0.6*(0.9-`X')^2                  
		if `k'==5 gen m5 = 0.2 - 0.6*(0.9-`X')^2 + sqrt(1-`X') - 0.1*`X'*cos(30*`X')
		if `k'==6 gen m6 = 0.4 + 0.25*sin(8*`X'-5) + 0.4*exp(-16*(4*`X'-2.5)^2)
		gen y`k' = 0*T + m`k' + e
		}
    // terms for estimating outcome equation
		gen lin  = `X'                                     
		gen sq   = (0.9-`X')^2 
		gen cube_p_tenth = -5*(`X'-0.7)^3 - 10*(`X'-0.6)^10 
		gen sqrt = sqrt(1-`X') 
		gen exp  = exp(-200*(`X'-0.7)^2)           
		gen exp2 = exp(-16*(4*`X'-2.5)^2)
		gen sqrt_p_cos = sqrt(1-`X') - 0.1*`X'*cos(30*`X')
		gen sin  = sin(8*`X'-5)
		gen pX2 = `X'^2
	// regression specs
	    global reg1 = "lin"
	    global reg2 = "lin exp"
	    global reg3 = "sq cube_p_tenth"
	    global reg4 = "sq sqrt"
	    global reg5 = "sq sqrt_p_cos"
	    global reg6 = "sin exp2"
		
** Organize
*drop  e u Tstar
aorder
}
end

capture program drop fouriercoeffs
program define fouriercoeffs
#delimit ;
matrix K = ( 3.8045830, -1.0764500, 0.2052452, 0.0000000, 0.0000000, 0.0357783, -0.0204930, 0.0052849, 0.0000000, 0.0000000,-0.0208564 \
	         1.6826000, -0.1599060, 0.0739300, 0.0000000, 0.0000000, 0.0077514, -0.0056401, 0.0016085, 0.0000000, 0.0000000,-0.0036519 \
	         0.8935404, -0.1201665, 0.1072360,-0.0283796, 0.0066526, 0.0561901, -0.0460486, 0.0295689, -.0151546, 0.0053763,-0.0291041 \
	         2.1070840, -0.4891177, 0.0894879, 0.0000000, 0.0000000, 0.6159659,  0.2563425,-0.0641854, 0.0000000, 0.0000000,-2.1959240 \
	         2.0815430, -0.4648486, 0.0810960, 0.0000000, 0.0000000,-0.5374329, -0.2979904, 0.0745325, 0.0000000, 0.0000000, 2.1522620 )';
#delimit cr   
end


*frolichlogitdgp, n(100) design(1) curves(1 2 3 4 5 6) var(0.1)
